R Code for Hoefer et al. 2024 - Sensors vs Surveyors: Comparing passive acoustic monitoring, camera trapping and observer-based monitoring for terrestrial mammals.

Author

Sebastian Hoefer

Published

October 30, 2024

Load packages

Load packages needed for data manipulation and analyses.

library(tidyverse)    # Data wrangling and visualization
library(colorspace)   # Color manipulation
library(cmdstanr)     # Interface to 'CmdStan'
library(brms)         # Bayesian regression models
library(rstan)        # Interface to Stan
library(ggeffects)    # Marginal effects plots
library(DHARMa)       # Residual diagnostics
library(emmeans)      # Estimated marginal means
library(tidybayes)    # Tidy data for Bayesian models
library(vegan)        # Community ecology analysis
library(EcolUtils)    # Ecology utilities
library(patchwork)    # Combine ggplots
library(HDInterval)   # HPD intervals
library(report)       # Automated reporting
library(ggsignif)     # Significance levels in ggplot
library(ggimage)      # Images in ggplot
library(cowplot)      # Plot composition in ggplot
library(paletteer)    # For colour palettes

Functions

Custom ggplot theme for visualisation.

my.theme <- function(){
  theme_classic() +
    theme(text = element_text(family = "Avenir Next"),
          axis.title.y = element_text(margin = margin(t = 0,r = 20,b = 0,l = 0)),
          axis.title.x = element_text(margin = margin(t = 20,r = 0,b = 0,l = 0)), 
          plot.margin = unit(c(5, 10, 5, 10), units = "mm"),
          strip.background = element_rect(fill = "#CCCCFF"),
          strip.text.x = element_text(size = 20),
          axis.title = element_text(size = 20),
          axis.text = element_text(size = 18),
          legend.text = element_text(size = 15),
          legend.title = element_text(size = 15))
}

Load data

All mammal data

mammal_data <- read.csv("mammal_data.csv")
mammal_data_vocal <- mammal_data %>% 
  filter(vocal == "yes",
         !common.name %in% c("Yellow-bellied Sheathtail-bat")) # Removed because mostly out of human audible range

Function to summarise the total number of species per method for each plot across all sites. Adding missing zero values for methods that didn’t capture any species at a given plot.

summarise_mammal_data <- function(df) {
  
  site_plots <- c("Tarcutta.DryA", "Tarcutta.DryB", "Tarcutta.WetA", "Tarcutta.WetB", 
                  "Duval.DryA", "Duval.DryB", "Duval.WetA", "Duval.WetB", 
                  "Mourachan.DryA", "Mourachan.DryB", "Mourachan.WetA", "Mourachan.WetB",
                  "Wambiana.DryA", "Wambiana.DryB", "Wambiana.WetA", "Wambiana.WetB",
                  "Undara.DryA", "Undara.DryB", "Undara.WetA", "Undara.WetB",
                  "Rinyirru.DryA", "Rinyirru.DryB", "Rinyirru.WetA", "Rinyirru.WetB")

  sites <- c("Tarcutta", "Duval", "Mourachan", "Wambiana", "Undara", "Rinyirru")
  
  methods <- c("PAM.all", "OBM", "PAM.survey", "camera")
  
  df <- df %>%
    unite(site.plot, c(site, plot), sep = ".", remove = FALSE) %>%
    select(-plot) %>%
    group_by(site, site.plot, assessment.method) %>%
    summarise(richness = length(unique(common.name))) %>% # Needs to be common.name because of Dingo and Dog having the same scientific.name
    ungroup() %>%
    select(site, site.plot, assessment.method, richness) %>%
    full_join(crossing(site.plot = unique(.$site.plot),
                       assessment.method = unique(.$assessment.method)),
              by = c("site.plot", "assessment.method")) %>%
    replace_na(list(richness = 0)) %>%
    mutate(site.plot = factor(site.plot, levels = site_plots),
           site = str_extract(site.plot, "^[^.]*"),
           site = factor(site, levels = sites),
           assessment.method = factor(assessment.method, levels = methods),
           richness = as.numeric(richness))

  return(df)
}

Run function.

mammal_data_summary <- summarise_mammal_data(mammal_data)
mammal_data_vocal_summary <- summarise_mammal_data(mammal_data_vocal)

PAM data

pam_data <- read_csv("mammals_PAM_daily_all.csv") %>% 
  filter(MANUAL.ID != "UNSURE",
         common.name != "Ghost Bat") %>% 
  mutate(MANUAL.ID = ifelse(MANUAL.ID == TRUE, 1, 0))
Rows: 198443 Columns: 26
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (15): IN FILE, INDIR, Recording, AudioLink, MANUAL.ID, template, site, plot, common.name, scientific.name, fami...
dbl   (8): CHANNEL, OFFSET, DURATION, Start, End, RecordingID, eucledian.distance, SamplingDay
dttm  (1): date.time
date  (1): date
time  (1): time

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

BirdNET Embeddings Performance Check

Custom function

Created function that creates a model for each template.

run_model <- function(df) {
  model_logistic <- glm(MANUAL.ID ~ eucledian.distance, family = binomial(link = "logit"), data = df)
  df$prediction <- predict(model_logistic, type = "response")
  return(df)
}

# Apply the model to each template and combine the results
all_data <- pam_data %>%
  group_by(template) %>%
  group_modify(~run_model(.x)) %>%
  ungroup()

Visualisation

# Create a single plot with facets for each template
(combined_plot <- ggplot(all_data, aes(x = eucledian.distance)) +  
  geom_point(aes(y = MANUAL.ID), alpha = 0.2) +
  geom_line(aes(y = prediction), col = "coral2") +
  facet_wrap(~ template, scales = "free_x") +
  labs(x = "Euclidean Distance", 
       y = "Probability of correct prediction") +
  scale_y_continuous(limits = c(0,1)) +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        axis.title = element_text(size = 18, face = "bold"),
        axis.title.x = element_text(margin = margin(t = 25, r = 0, b = 0, l = 0)),
        axis.title.y = element_text(margin = margin(t = 0, r = 25, b = 0, l = 0)),
        plot.margin = unit(c(5, 10, 5, 10), units = "mm")))

Species Richness

Given the inherent inability of passive acoustic monitoring to detect non-vocalising species, the mammal community was considered in two ways: i) entire mammal community (All Mammals), and ii) vocal mammals only (Vocal Mammals).

All Mammals

Considering the entire mammal community with vocal and non-vocal species.

Bayesian Analysis

Priors

Defining priors for Bayesian models.

Fit model

Compare Models

Diagnostics

Investigation

Methods pairwise comparison across sites
(diff.methods.avg <- methods.brm1b %>%
  emmeans(~assessment.method) %>%
  regrid(transform = "none") %>%
  pairs() %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n())) #to see if there is any change
# A tibble: 6 × 5
  contrast             `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
  <chr>                                      <dbl>       <dbl>       <dbl>                       <dbl>
1 OBM - PAM.survey                           2.46        1.85        3.12                       1     
2 OBM - camera                               1.44        1.16        1.77                       0.999 
3 PAM.all - OBM                              0.863       0.691       1.04                       0.0729
4 PAM.all - PAM.survey                       2.13        1.62        2.72                       1     
5 PAM.all - camera                           1.25        0.992       1.52                       0.977 
6 PAM.survey - camera                        0.587       0.430       0.744                      0     
(diff.methods.avg.rev <- methods.brm1b %>%
  emmeans(~assessment.method) %>%
  regrid(transform = "none") %>%
  pairs(reverse = TRUE) %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n())) #to see if there is any change
# A tibble: 6 × 5
  contrast             `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
  <chr>                                      <dbl>       <dbl>       <dbl>                       <dbl>
1 OBM - PAM.all                              1.16        0.932       1.39                     0.927   
2 PAM.survey - OBM                           0.406       0.304       0.513                    0       
3 PAM.survey - PAM.all                       0.469       0.359       0.601                    0       
4 camera - OBM                               0.693       0.555       0.848                    0.000833
5 camera - PAM.all                           0.801       0.651       0.999                    0.0229  
6 camera - PAM.survey                        1.70        1.26        2.17                     1       

Visualisation

pal <- c("#189F9F","#FF8C00", "#18709F", "#A034F0")
pal_dark <- after_scale(darken(pal, .1, space = "HLS"))

(methods.boxplot <- ggplot(mammal_data_summary, aes(assessment.method, richness)) +
  geom_boxplot(aes(color = assessment.method,
                   fill = after_scale(desaturate(lighten(color, .8), .4))),
               position = position_dodge2(width = 0.6),
               width = .6, outlier.shape = NA) +
  geom_point(aes(fill = assessment.method), color = "transparent",
             position = position_dodge2(width = 0.6),
             shape = 21, stroke = .4, size = 3.5, alpha = .3) +
  my.theme() +
  scale_y_continuous(limits = c(0,20),breaks = seq(0, 20, by = 2),
                     name = "Total Species Richness") +
  scale_x_discrete(name = "", labels = c("", "", "", "")) +
  scale_fill_manual(values = pal, guide = "none") +
  scale_colour_manual(values = pal_dark, name = "", labels = c("PAM (all audio)","OBM", "PAM (survey period)","Camera Trap")) +
  theme(legend.text = element_text(size = 14, margin = margin(t = 5)),
        legend.position = c(y=0.45, x=-0.15),
        legend.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(5, 10, 20, 10), units = "mm")) +
  guides(colour = guide_legend(nrow = 1)) +
    geom_signif(comparisons = list(c("PAM.all", "PAM.survey"),
                                   c("OBM", "PAM.survey"),
                                   c("OBM", "camera"),
                                   c("PAM.survey", "camera")),
                annotations = c("2.13, HDI 95%[1.62,2.72]",
                                "2.46, HDI 95%[1.85,3.12]",
                                "1.44, HDI 95%[1.16,1.77]",
                                "0.59, HDI 95%[0.43,0.74]"),
                y_position = c(16, 17.5, 19, 12)))

Vocal Mammals

Considering vocal mammals only and excluding non-vocal mammals from the analysis.

Bayesian Analysis

Priors

Defining priors for Bayesian models.

Fit model

Compare Models
Diagnostics

Investigation

Methods pairwise comparison across sites
(diff.methods.vocal.avg <- methods.vocal.brm1b %>%
  emmeans(~assessment.method) %>%
  regrid(transform = "none") %>%
  pairs() %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n())) #to see if there is any change
# A tibble: 6 × 5
  contrast             `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
  <chr>                                      <dbl>       <dbl>       <dbl>                       <dbl>
1 OBM - PAM.survey                            1.34       0.989        1.74                       0.980
2 OBM - camera                                1.46       1.05         1.87                       0.995
3 PAM.all - OBM                               1.58       1.24         1.98                       1    
4 PAM.all - PAM.survey                        2.12       1.61         2.69                       1    
5 PAM.all - camera                            2.29       1.73         2.93                       1    
6 PAM.survey - camera                         1.08       0.785        1.43                       0.697
(diff.methods.vocal.avg.rev <- methods.vocal.brm1b %>%
  emmeans(~assessment.method) %>%
  regrid(transform = "none") %>%
  pairs(reverse = TRUE) %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n())) #to see if there is any change
# A tibble: 6 × 5
  contrast             `Average fractional change` `Lower HDI` `Upper HDI` `Probability of difference`
  <chr>                                      <dbl>       <dbl>       <dbl>                       <dbl>
1 OBM - PAM.all                              0.633       0.498       0.799                      0     
2 PAM.survey - OBM                           0.746       0.555       0.974                      0.0204
3 PAM.survey - PAM.all                       0.472       0.357       0.600                      0     
4 camera - OBM                               0.686       0.514       0.905                      0.005 
5 camera - PAM.all                           0.436       0.325       0.550                      0     
6 camera - PAM.survey                        0.923       0.681       1.23                       0.303 

Visualisation

pal <- c("#189F9F","#FF8C00", "#18709F","#A034F0")
pal_dark <- after_scale(darken(pal, .1, space = "HLS"))

(methods.boxplot.vocal <- ggplot(mammal_data_vocal_summary, aes(assessment.method, richness)) +
  geom_boxplot(aes(color = assessment.method,
                   fill = after_scale(desaturate(lighten(color, .8), .4))),
               position = position_dodge2(width = 0.6),
               width = .6, outlier.shape = NA) +
  geom_point(aes(fill = assessment.method), color = "transparent",
             position = position_dodge2(width = 0.6),
             shape = 21, stroke = .4, size = 3.5, alpha = .3) +
  my.theme() +
  scale_y_continuous(limits = c(0,20),breaks = seq(0, 20, by = 2),
                     name = "") +
  scale_x_discrete(name = "", labels = c("", "", "", "")) +
  scale_fill_manual(values = pal, guide = "none") +
  scale_colour_manual(values = pal_dark, name = "", labels = c("PAM (all audio)","OBM", "PAM (survey period)","Camera Trap")) +
  theme(legend.text = element_text(size = 14, margin = margin(t = 5)),
        legend.position = c(y=0.45, x=-0.15),
        legend.background = element_rect(fill = "transparent"),
        plot.margin = unit(c(5, 10, 20, 10), units = "mm")) +
  guides(colour = guide_legend(nrow = 1)) +
  geom_signif(comparisons = list(c("PAM.all", "OBM"),
                                 c("PAM.all", "PAM.survey"),
                                 c("PAM.all", "camera"),
                                 c("OBM", "camera")),
                annotations = c("1.58, HDI 95%[1.24,1.98]",
                                "2.12, HDI 95%[1.61,2.69]",
                                "2.29, HDI 95%[1.73,2.93]",
                                "1.46, HDI 95%[1.05,1.87]"),
                y_position = c(12.5, 14, 15.5, 9)))

Combined Figure
(methods.vocal.combined <- (methods.boxplot + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (methods.boxplot.vocal + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   theme(legend.position = c(-0.25,-0.2), legend.text = element_text(size = 22),
         plot.margin = unit(c(5, 10, 30, 10), units = "mm")))

Correlation PAM vocal mammals & OBM mammals

Investigate a correlation between species richness of vocal mammals detected by PAM and species richness of the entire mammal community detected by OBM.

cor_mammal_data <- mammal_data_summary |> 
  filter(assessment.method %in% c("OBM", "PAM.all")) |> 
  group_by(site.plot) |> 
  pivot_wider(names_from = "assessment.method", #transformed dataset from long to wide format
              values_from = "richness",
              values_fill = list(number=0))

(corr_plot <- ggplot(cor_mammal_data, aes(PAM.all,OBM)) +
  geom_point(aes(col = site), size = 3) +
  geom_smooth(method = "lm") +
  my.theme() +
  theme(legend.position = "bottom") +
  scale_y_continuous(name = "Species Richness OBM", breaks = seq(4, 16, by = 4)) +
  scale_x_continuous(name = "Species Richness PAM", breaks = seq(2, 13, by = 2)) +
  scale_color_manual(name = "", values = c("seagreen", "steelblue","salmon2", "darkorchid1", "aquamarine3", "goldenrod1")) +
  guides(color = guide_legend(nrow = 1)) +
  annotate("text", x = 5, y = 15,
           label = "R^2 == 0.21", parse = TRUE, size = 5))
`geom_smooth()` using formula = 'y ~ x'

# Perform Pearson's correlation test
correlation_test <- cor.test(cor_mammal_data$PAM.all, cor_mammal_data$OBM, method = "pearson")

print(correlation_test)

    Pearson's product-moment correlation

data:  cor_mammal_data$PAM.all and cor_mammal_data$OBM
t = 2.4133, df = 22, p-value = 0.02458
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.06636028 0.72677443
sample estimates:
      cor 
0.4575097 
(R.squared <- 0.4575097^2)
[1] 0.2093151

Community Composition

All Mammals

Considering the entire mammal community with vocal and non-vocal species.

Prepare data

mammal_community3 <- mammal_data %>% 
  filter(!recapture %in% c("y")) %>%  # Removed recaptures
  select(site, assessment.method,common.name) %>% # Selected relevant grouping variables
  group_by(across(everything())) %>% # Grouped by all available columns
  mutate(number = n()) %>% # Created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% # Merged the duplicate rows
  pivot_wider(names_from = "common.name", # Transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% # Replaced NAs with 0
  arrange(site) %>% # Ordered by site
  mutate(assessment.method = factor(assessment.method)) # Turned asssessment methods into factors

In preparation for community analyses, only keep columns and rows with values above 0.

all <- mammal_community3 %>% select(where(~ any(. != 0))) %>% filter(rowSums(.[,3:ncol(.)]) > 0)

Transformation of the data to binary (presence/absence) for Jaccard dissimilarity.

all_jacc <- all %>%
  select(where(is.numeric)) %>% # Removed columns with non-numerical values
  decostand(method="pa", MARGIN = 1) %>% # Standardised presence/absence (method = "pa") by rows (MARGIN = 1)
  cbind(all[,1:2]) %>% 
  select(site, assessment.method, everything())

NMDS Ordination

Jaccard dissimilarity

set.seed(11)
all_nmds_jacc <-  metaMDS(all_jacc[4:ncol(all_jacc)], distance="jaccard", k=3,trymax=100, autotransform = FALSE)
Run 0 stress 0.1055242 
Run 1 stress 0.1207629 
Run 2 stress 0.1047353 
... New best solution
... Procrustes: rmse 0.0249656  max resid 0.08737959 
Run 3 stress 0.105525 
Run 4 stress 0.1083786 
Run 5 stress 0.109782 
Run 6 stress 0.104735 
... New best solution
... Procrustes: rmse 0.0002099691  max resid 0.0007380798 
... Similar to previous best
Run 7 stress 0.1047354 
... Procrustes: rmse 0.0002734112  max resid 0.0009616091 
... Similar to previous best
Run 8 stress 0.114386 
Run 9 stress 0.1207632 
Run 10 stress 0.1133275 
Run 11 stress 0.1097819 
Run 12 stress 0.1123875 
Run 13 stress 0.1133275 
Run 14 stress 0.1145607 
Run 15 stress 0.1133275 
Run 16 stress 0.1055252 
Run 17 stress 0.1166412 
Run 18 stress 0.1145605 
Run 19 stress 0.1155145 
Run 20 stress 0.1200923 
*** Best solution repeated 2 times

Visualisation

Jaccard dissimilarity

Function to create nMDS plots with centroids.

create_plot <- function(data, nmds_data) {
  methods <- nmds_data$points
  species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
  df <- cbind.data.frame(data[,1:2], methods)
  gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
    mutate(assessment.method = factor(assessment.method, levels = c("PAM.all", "OBM","PAM.survey", "camera")))
  ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
    geom_point(size=3) + # Smaller outside points
    geom_point(aes(x=mean.x,y=mean.y),size=5) + # This controls the centroids
    scale_colour_manual(values = c("#189F9F","#FF8C00", "#18709F","#A034F0"),
                        name = "",
                        labels = c("PAM (all audio)", "OBM", "PAM (survey period)","Camera Trap")) +
    geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) + # This controls the connections between points and centroids
    my.theme() +
    theme(legend.position = "bottom")
}

(nmds.all.plot.jacc.centroid <- create_plot(all_jacc, all_nmds_jacc))

Investigation

Adonis

Jaccard dissimilarity
all.dist2 <- vegdist(all_jacc[3:ncol(all_jacc)], 'jaccard')
(all.adonis <- adonis2(all.dist2 ~ assessment.method, data=all_jacc))
Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 999

adonis2(formula = all.dist2 ~ assessment.method, data = all_jacc)
         Df SumOfSqs      R2     F Pr(>F)    
Model     3   1.9870 0.30827 2.971  0.001 ***
Residual 20   4.4586 0.69173                 
Total    23   6.4456 1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Permutation test: significant difference between assessment methods for mammal communities.

#Let's which method is different from which
(adonis.methods.jacc <- adonis.pair(all.dist2, all_jacc$assessment.method, corr.method = "holm"))
             combination SumsOfSqs   MeanSqs  F.Model        R2     P.value P.value.corrected
1         camera <-> OBM 0.4559165 0.4559165 1.692142 0.1447247 0.105894106       0.211788212
2     camera <-> PAM.all 0.9499132 0.9499132 4.958432 0.3314807 0.000999001       0.005994006
3  camera <-> PAM.survey 1.1114297 1.1114297 4.864324 0.3272483 0.003996004       0.015984016
4        OBM <-> PAM.all 0.5417048 0.5417048 2.492023 0.1994891 0.018981019       0.056943057
5     OBM <-> PAM.survey 0.6990302 0.6990302 2.748990 0.2156241 0.002997003       0.014985015
6 PAM.all <-> PAM.survey 0.2159058 0.2159058 1.223748 0.1090320 0.290709291       0.290709291
#comparing all methods with each other

Vocal Mammals

Considering vocal mammals only and excluding non-vocal mammals from the analysis.

Prepare data

mammal_community4 <- mammal_data_vocal %>% 
  filter(!recapture %in% c("y")) %>%  
  select(site, assessment.method,common.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>% 
  mutate(assessment.method = factor(assessment.method))
all_vocal <- mammal_community4 %>% select(where(~ any(. != 0))) %>% filter(rowSums(.[,3:ncol(.)]) > 0)
#Function decostand() standardises presence/absence (method = "pa") by rows (MARGIN = 1)
all_vocal_jacc <- all_vocal %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1) %>% 
  cbind(all_vocal[,1:2]) %>% 
  select(site, assessment.method, everything())

NMDS Ordination

Jaccard dissimilarity

set.seed(11)
all_nmds_jacc_vocal <-  metaMDS(all_vocal_jacc[4:ncol(all_vocal_jacc)], distance="jaccard", k=3,trymax=100, autotransform = FALSE)

Visualisation

Jaccard dissimilarity

create_plot <- function(data, nmds_data) {
  methods <- nmds_data$points
  species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
  df <- cbind.data.frame(data[,1:2], methods)
  gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
    mutate(assessment.method = factor(assessment.method, levels = c("PAM.all", "OBM","PAM.survey", "camera")))
  ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
    geom_point(size=3) +
    geom_point(aes(x=mean.x,y=mean.y),size=5) +
    scale_colour_manual(values = c("#189F9F","#FF8C00", "#18709F","#A034F0"),
                        name = "",
                        labels = c("PAM (all audio)", "OBM", "PAM (survey period)","Camera Trap")) +
    geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) +
    my.theme() + theme(legend.position = "bottom")
}

(nmds.all.vocal.plot.jacc.centroid <- create_plot(all_vocal_jacc, all_nmds_jacc_vocal))

Investigation

Adonis

Jaccard dissimilarity
all.dist_vocal <- vegdist(all_vocal_jacc[3:ncol(all_vocal_jacc)], 'jaccard')
(all.adonis_vocal <- adonis2(all.dist_vocal ~ assessment.method, data=all_vocal_jacc))
Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 999

adonis2(formula = all.dist_vocal ~ assessment.method, data = all_vocal_jacc)
         Df SumOfSqs     R2      F Pr(>F)   
Model     3   1.3414 0.2559 2.2927   0.01 **
Residual 20   3.9006 0.7441                 
Total    23   5.2420 1.0000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Permutation test: significant difference between assessment methods for mammal communities.

#Let's which method is different from which
(adonis.methods.jacc_vocal <- adonis.pair(all.dist_vocal, all_vocal_jacc$assessment.method, corr.method = "holm"))
             combination SumsOfSqs   MeanSqs   F.Model         R2     P.value P.value.corrected
1         camera <-> OBM 0.4621827 0.4621827 2.1635228 0.17786976 0.071928072        0.28771229
2     camera <-> PAM.all 0.6947730 0.6947730 4.1757138 0.29456815 0.004995005        0.02497502
3  camera <-> PAM.survey 0.8980116 0.8980116 4.4172857 0.30638816 0.003996004        0.02397602
4        OBM <-> PAM.all 0.1200808 0.1200808 0.6429681 0.06041248 0.713286713        0.71328671
5     OBM <-> PAM.survey 0.2919053 0.2919053 1.3050666 0.11544086 0.233766234        0.70129870
6 PAM.all <-> PAM.survey 0.2159058 0.2159058 1.2237477 0.10903200 0.290709291        0.70129870
#comparing all methods with each other

Combined Figure

(nmds.jacc.combined <- (nmds.all.plot.jacc.centroid + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (nmds.all.vocal.plot.jacc.centroid + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   theme(legend.position = c(-0.25,-0.2), legend.text = element_text(size = 22),
         plot.margin = unit(c(5, 10, 30, 10), units = "mm")) +
   guides(colour = guide_legend(nrow = 1)))

Survey Effort (Species accumulation curves)

All Mammals

Considering the entire mammal community with vocal and non-vocal species.

Prepare data

Merged duplicate rows for duplicate observations and transformed the dataset to a wide format.

mammal_community <- mammal_data %>% 
  select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>% 
  filter(assessment.method != "PAM.all") %>% 
  mutate(SamplingDay = as.factor(SamplingDay))

Added a category for the total richness.

result.i <- vector("list",length(unique(mammal_community$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(mammal_community$site))){ #looped over each study site
  site.i <- mammal_community[mammal_community$site == unique(mammal_community$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
  for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
    
    
    SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] #subset to the SamplingDay
    
    meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("total"))
    
    SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
    
    colnames(total.a) <- colnames(mammal_community)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.i <- do.call("rbind.data.frame",result.i)

mammal_community2 <- rbind.data.frame(mammal_community, result.i)

Added a remote sensing category to assessment methods. Remote sensing was a combination of PAM and camera trapping.

result.j <- vector("list",length(unique(mammal_community$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(mammal_community$site))){ #looped over each study site
  site.i <- mammal_community[mammal_community$site == unique(mammal_community$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
  for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
    
    
    SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] %>% 
      filter(assessment.method %in% c("camera", "PAM.survey"))
    
    meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("remote"))
    
    SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
    
    colnames(total.a) <- colnames(mammal_community)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.j[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.j <- do.call("rbind.data.frame",result.j)

mammal_community2 <- rbind.data.frame(mammal_community2, result.j) %>% 
  drop_na()

Loop

Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.

# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 30, ncol = (ncol(mammal_community2)))) %>% 
  replace(is.na(.), 0) %>% 
  set_names(names(mammal_community2)) %>% 
  select(-(1:3))

result.i <- vector("list",length(unique(mammal_community2$site))) #created an empty list that was filled as the loop ran

set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(mammal_community2$site))){ #looped over each study site
  site.i <- mammal_community2[mammal_community2$site == unique(mammal_community2$site)[i],] #subset to the site
  
  dates.i <- ifelse(unique(site.i$site) %in% c("Tarcutta", "Undara", "Wambiana", "Rinyirru"), 28, 21) #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day 
  
  result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
  for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
    
    
    method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
    
    if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
      sampling.a <- method.a$assessment.method[1] #saved the name of the method
      
      method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
      
      if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
        method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
      
      accum.a <- specaccum(method.a,"random") #calculated data
      
      #created data frame of metadata, richness, and standard deviation
      accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
      colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
      
      accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
      accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
      
      result.a[[a]] <- accum.a}} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

Transformed resulting list to a data frame.

mammals_wide.result <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data frame
mutate(assessment.method = factor(assessment.method, levels = c("OBM",
                                                                "camera",
                                                                "PAM.survey",
                                                                "remote",
                                                                "total")),
         site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

Visualisation

Plot of rarefaction curves for each method split by sites.

(specaccum_mammals <- ggplot(mammals_wide.result, aes(x = day, y = richness, col = assessment.method)) + 
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.3))),
              linetype = 0.1) +
  geom_vline(xintercept=c(7,14,21, 28), linetype="dotted", colour = "black") +
  geom_line(aes(group = assessment.method),linewidth = 1.3) +
  my.theme() +
  facet_wrap(~site) +
  theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
        strip.background = element_rect(fill = "lightgrey"),
        axis.title = element_text(size = 26),
        axis.text = element_text(size = 26),
        legend.text = element_text(size = 26),
        legend.position = "bottom",
        strip.text.x = element_text(size = 26)) +
  scale_y_continuous(limits = c(0,25), breaks = seq(0, 25, by = 5),
                     name = "Species Richness") +
  scale_x_continuous(limits = c(0,28), breaks = seq(0, 30, by = 7),
                     name = "Survey Effort (Days)") +
  scale_colour_manual(values = c("#FF8C00", "#A034F0", "#189F9F","#D14103","black"),
                      name = "",
                      labels = c("OBM","Camera Trap", "PAM (survey period)", "Remote Sensing",
                                 "Total Richness")))

Vocal Mammals

Considering vocal mammals only and excluding non-vocal mammals from the analysis.

Prepare data

Merged duplicate rows for duplicate observations and transformed the dataset to a wide format.

mammal_community_vocal <- mammal_data_vocal %>% 
  select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>% 
  filter(assessment.method != "PAM.all") %>% 
  mutate(SamplingDay = as.factor(SamplingDay))

Added a total category to assessment methods.

result.i <- vector("list",length(unique(mammal_community_vocal$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(mammal_community_vocal$site))){ #looped over each study site
  site.i <- mammal_community_vocal[mammal_community_vocal$site == unique(mammal_community_vocal$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
  for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
    
    
    SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] #subset to the SamplingDay
    
    meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("total"))
    
    SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
    
    colnames(total.a) <- colnames(mammal_community_vocal)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.i <- do.call("rbind.data.frame",result.i)

mammal_community_vocal2 <- rbind.data.frame(mammal_community_vocal, result.i)

Added a remote sensing category to assessment methods. Remote sensing was a combination of PAM and camera trapping.

result.j <- vector("list",length(unique(mammal_community_vocal$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(mammal_community_vocal$site))){ #looped over each study site
  site.i <- mammal_community_vocal[mammal_community_vocal$site == unique(mammal_community_vocal$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$SamplingDay))) #created an empty list to be filled by each SamplingDay
  for(a in 1:length(unique(site.i$SamplingDay))){ #looped over each SamplingDay
    
    
    SamplingDay.a <- site.i[site.i$SamplingDay == unique(site.i$SamplingDay)[a],] %>% 
      filter(assessment.method %in% c("camera", "PAM.survey"))
    
    meta.a <- cbind.data.frame(SamplingDay.a$site[1],SamplingDay.a$SamplingDay[1],c("remote"))
    
    SamplingDay.a <- SamplingDay.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(SamplingDay.a))) #calculated column sums
    
    colnames(total.a) <- colnames(mammal_community_vocal)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.j[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.j <- do.call("rbind.data.frame",result.j)

mammal_community_vocal2 <- rbind.data.frame(mammal_community_vocal2, result.j) %>% 
  drop_na()

Loop

Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.

# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 30, ncol = (ncol(mammal_community_vocal2)))) %>% 
  replace(is.na(.), 0) %>% 
  set_names(names(mammal_community_vocal2)) %>% 
  select(-(1:3))

result.i <- vector("list",length(unique(mammal_community_vocal2$site))) #created an empty list that was filled as the loop ran

set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(mammal_community_vocal2$site))){ #looped over each study site
  site.i <- mammal_community_vocal2[mammal_community_vocal2$site == unique(mammal_community_vocal2$site)[i],] #subset to the site
  
  dates.i <- ifelse(unique(site.i$site) %in% c("Tarcutta", "Undara", "Wambiana", "Rinyirru"), 28, 21) #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day 
  
  result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
  for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
    
    
    method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
    
    if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
      sampling.a <- method.a$assessment.method[1] #saved the name of the method
      
      method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
      
      if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
        method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
      
      accum.a <- specaccum(method.a,"random") #calculated data
      
      #created data frame of metadata, richness, and standard deviation
      accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
      colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
      
      accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
      accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
      
      result.a[[a]] <- accum.a}} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

Transformed resulting list to a data frame.

mammals_wide.result_vocal <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data frame
mutate(assessment.method = factor(assessment.method, levels = c("OBM",
                                                                "camera",
                                                                "PAM.survey",
                                                                "remote",
                                                                "total")),
         site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

Visualisation

Plot of rarefaction curves for each method split by sites.

(specaccum_mammals_vocal <- ggplot(mammals_wide.result_vocal, aes(x = day, y = richness, col = assessment.method)) + 
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.3))),
              linetype = 0.1) +
  geom_vline(xintercept=c(7,14,21, 28), linetype="dotted", colour = "black") +
  geom_line(aes(group = assessment.method), linewidth = 1.3) +
  my.theme() +
  facet_wrap(~site) +
  theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
        strip.background = element_rect(fill = "lightgrey"),
        axis.title = element_text(size = 26),
        axis.text = element_text(size = 26),
        legend.text = element_text(size = 26),
        legend.position = "bottom",
        strip.text.x = element_text(size = 26)) +
  scale_y_continuous(limits = c(0,15), breaks = seq(0, 15, by = 5),
                     name = "") +
  scale_x_continuous(limits = c(0,28), breaks = seq(0, 30, by = 7),
                     name = "Survey Effort (Days)") +
  scale_colour_manual(values = c("#FF8C00", "#A034F0", "#189F9F","#D14103" ,"black"),
                      name = "",
                      labels = c("OBM","Camera Trap", "PAM (survey period)", "Remote Sensing",
                                 "Total Richness")))

Combined Figure

(specaccum_mammals_combined <- (specaccum_mammals + ggtitle('All mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (specaccum_mammals_vocal + ggtitle('Vocal mammals') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   theme(legend.position = c(-0.1,-0.2), legend.text = element_text(size = 22),
         plot.margin = unit(c(5, 10, 30, 10), units = "mm")) +
   guides(colour = guide_legend(nrow = 1)))

PAM survey length

How many days until PAM.all reaches values similar to OBM for all mammals and only vocal mammals. In other words how long should we leave recorders out on average until we get the best number of species.

all_mammals <- mammal_data %>% 
  filter(assessment.method == "PAM.all") %>% 
  select(site, SamplingDay, assessment.method,common.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "common.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>%
  mutate(SamplingDay = as.factor(SamplingDay))

Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.

# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 1235, ncol = (ncol(all_mammals)))) %>% 
  replace(is.na(.), 0) %>% 
  set_names(names(all_mammals)) %>% 
  select(-(1:3))

result.i <- vector("list",length(unique(all_mammals$site))) #created an empty list that was filled as the loop ran

set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(all_mammals$site))){ #looped over each study site
  site.i <- all_mammals[all_mammals$site == unique(all_mammals$site)[i],] #subset to the site
  
   #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day 
  site_values <- c(Tarcutta = 986, Duval = 679, Mourachan = 783, Wambiana = 1235, Undara = 780, Rinyirru = 488)
  
  dates.i <- site_values[match(unique(site.i$site), names(site_values))]
  
  result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
  for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
    
    
    method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
    
    if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
      sampling.a <- method.a$assessment.method[1] #saved the name of the method
      
      method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
      
      if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
        method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
      
      accum.a <- specaccum(method.a,"random") #calculated data
      
      #created data frame of metadata, richness, and standard deviation
      accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
      colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
      
      accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
      accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
      
      result.a[[a]] <- accum.a}} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

Transformed resulting list to a data frame.

all_mammals.result <- do.call("rbind.data.frame",result.i) %>%  #compressed result.i list into data
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))
pam.max <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == max(day)) %>% 
  ungroup() %>% 
  select(site, richness.max = richness)

pam.28 <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == 28) %>% 
  ungroup() %>% 
  select(richness.28 = richness) %>% 
  mutate(prop.28 = round(richness.28/pam.max$richness.max, digits = 2))

pam.90 <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == 90) %>% 
  ungroup() %>% 
  select(richness.90 = richness) %>% 
  mutate(prop.90 = round(richness.90/pam.max$richness.max, digits = 2))

pam.180 <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == 180) %>% 
  ungroup() %>% 
  select(richness.180 = richness) %>% 
  mutate(prop.180 = round(richness.180/pam.max$richness.max, digits = 2))

pam.365 <- all_mammals.result %>% 
  group_by(site) %>% 
  filter(assessment.method == "PAM.all",
         day == 365) %>% 
  ungroup() %>% 
  select(richness.365 = richness) %>% 
  mutate(prop.365 = round(richness.365/pam.max$richness.max, digits = 2))

pam.rich.prop <- cbind(pam.max,pam.28,pam.90,pam.180,pam.365) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

(sampling.table <- pam.rich.prop %>% 
  select(Site = site, "PAM (total richness)" = richness.max,
         "PAM (proportion 28 days)" = prop.28,
         "PAM (proportion 90 days)" = prop.90,
         "PAM (proportion 180 days)" = prop.180,
         "PAM (proportion 365 days)" = prop.365) %>% 
  add_row(Site = "Average", "PAM (total richness)" = round(mean(pam.rich.prop$richness.max)),
          "PAM (proportion 28 days)" = round(mean(pam.rich.prop$prop.28), digits = 2),
          "PAM (proportion 90 days)" = round(mean(pam.rich.prop$prop.90), digits = 2),
          "PAM (proportion 180 days)" = round(mean(pam.rich.prop$prop.180), digits = 2),
          "PAM (proportion 365 days)" = round(mean(pam.rich.prop$prop.365), digits = 2)))
       Site PAM (total richness) PAM (proportion 28 days) PAM (proportion 90 days) PAM (proportion 180 days)
1     Duval                   10                     0.72                     0.86                      0.91
2 Mourachan                    7                     0.68                     0.86                      0.89
3  Rinyirru                    9                     0.61                     0.79                      0.87
4  Tarcutta                   13                     0.66                     0.78                      0.86
5    Undara                    7                     0.50                     0.75                      0.86
6  Wambiana                   10                     0.79                     0.91                      0.96
7   Average                    9                     0.66                     0.82                      0.89
  PAM (proportion 365 days)
1                      0.95
2                      0.92
3                      0.96
4                      0.94
5                      0.92
6                      0.99
7                      0.95

Visualisation

vline_data <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
                         xintercept = c(12, 175, 17, 88, 676,320)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

vline_data2 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
                         xintercept = c(12, 35, 12, 17, 12,9)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

hline_data_vocal <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(4.1, 6, 7, 6,12,10)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

hline_data_vocal2 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(3.9, 4, 6, 4,6,6)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))



hline_data_cam <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(4, 4, 6, 4,6,6)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

vline_data_cam <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
                         xintercept = c(12, 35, 12, 17, 12,9)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

hline_data_28 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(5.48, 3.51, 7.88, 4.32,8.52,6.87)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

vline_data_28 <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Duval", "Tarcutta"),
                         xintercept = c(28, 28, 28, 21, 21,28)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

hline_data_total <- data.frame(site = c("Rinyirru", "Undara", "Wambiana", "Mourachan", "Tarcutta", "Duval"),
                               yintercept = c(8, 7, 11, 8,13,11)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

vline_data_total <- data.frame(site = c("Rinyirru", "Undara","Tarcutta"),
                         xintercept = c(210, 771, 968)) %>% 
  mutate(site = factor(site, levels = c("Rinyirru", "Undara","Tarcutta")))


(spec_all_mammals_plot <- ggplot(all_mammals.result, aes(x = day, y = richness, col = assessment.method)) +
    geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.5))),
              linetype = 0.1) +
    geom_line(aes(group = assessment.method), linewidth = 1.5) +
    my.theme() +
    facet_wrap(~site, scales = "free_x") +
    geom_hline(data = hline_data_vocal, aes(yintercept = yintercept), linetype="dashed", color = "#FF8C00",
               linewidth = 1) +
    geom_hline(data = hline_data_vocal2, aes(yintercept = yintercept), linetype="dashed", color = "#A034F0",
               linewidth = 1) +
    #geom_vline(data = vline_data, aes(xintercept = xintercept), linetype="dotted", color = "#FF8C00") +
    #geom_vline(data = vline_data2, aes(xintercept = xintercept), linetype="dotted", color = "#A034F0") +
    geom_text(data = vline_data, aes(x = 100, y = 14, label = paste0("DTM: ",xintercept)),
                                     vjust = -1, size = 8,inherit.aes = FALSE, color = "#FF8C00") +
    geom_text(data = vline_data2, aes(x = 100, y = 13, label = paste0("DTM: ",xintercept)),
                                     vjust = -1, size = 8,inherit.aes = FALSE, color = "#A034F0") +
    theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
        strip.background = element_rect(fill = "lightgrey"),
        axis.title = element_text(size = 26),
        axis.text = element_text(size = 26),
        legend.text = element_text(size = 26),
        legend.position = "bottom",
        strip.text.x = element_text(size = 26)) + 
    scale_y_continuous(limits = c(0,15), breaks = seq(0, 15, by = 5),
                     name = "Species Richness") +
    scale_x_continuous(name = "Survey Effort (Days)") +
    scale_colour_manual(values = c("#189F9F"),
                      name = "",
                      labels = c("PAM (all audio)")))

Session Info

sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS 15.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Brisbane
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] paletteer_1.6.0     cowplot_1.1.3       ggimage_0.3.3       ggsignif_0.6.4      report_0.5.9       
 [6] HDInterval_0.2.4    patchwork_1.3.0     EcolUtils_0.1       vegan_2.6-8         lattice_0.22-6     
[11] permute_0.9-7       tidybayes_3.0.7     emmeans_1.10.4      DHARMa_0.4.6        ggeffects_1.7.1    
[16] rstan_2.32.6        StanHeaders_2.32.10 brms_2.21.0         Rcpp_1.0.13         cmdstanr_0.8.0.9000
[21] colorspace_2.1-1    lubridate_1.9.3     forcats_1.0.0       stringr_1.5.1       dplyr_1.1.4        
[26] purrr_1.0.2         readr_2.1.5         tidyr_1.3.1         tibble_3.2.1        ggplot2_3.5.1      
[31] tidyverse_2.0.0    

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1     rstudioapi_0.16.0    jsonlite_1.8.9       datawizard_0.12.3    magrittr_2.0.3      
  [6] estimability_1.5.1   magick_2.8.5         farver_2.1.2         nloptr_2.1.1         rmarkdown_2.28      
 [11] fs_1.6.4             vctrs_0.6.5          minqa_1.2.8          htmltools_0.5.8.1    distributional_0.5.0
 [16] haven_2.5.4          gridGraphics_0.5-1   htmlwidgets_1.6.4    plyr_1.8.9           TMB_1.9.15          
 [21] mime_0.12            iterators_1.0.14     lifecycle_1.0.4      pkgconfig_2.0.3      gap_1.6             
 [26] Matrix_1.7-0         R6_2.5.1             fastmap_1.2.0        shiny_1.9.1          rbibutils_2.2.16    
 [31] digest_0.6.37        numDeriv_2016.8-1.1  rematch2_2.1.2       ps_1.8.0             qgam_1.3.4          
 [36] labeling_0.4.3       fansi_1.0.6          timechange_0.3.0     abind_1.4-8          mgcv_1.9-1          
 [41] compiler_4.4.1       doParallel_1.0.17    bit64_4.0.5          withr_3.0.1          backports_1.5.0     
 [46] inline_0.3.19        QuickJSR_1.3.1       pkgbuild_1.4.4       MASS_7.3-60.2        loo_2.8.0           
 [51] tools_4.4.1          httpuv_1.6.15        glue_1.7.0           promises_1.3.0       nlme_3.1-164        
 [56] grid_4.4.1           checkmate_2.3.2      reshape2_1.4.4       cluster_2.1.6        generics_0.1.3      
 [61] glmmTMB_1.1.9        gtable_0.3.5         tzdb_0.4.0           data.table_1.16.0    hms_1.1.3           
 [66] utf8_1.2.4           foreach_1.5.2        pillar_1.9.0         ggdist_3.3.2         yulab.utils_0.1.7   
 [71] vroom_1.6.5          later_1.3.2          posterior_1.6.0      splines_4.4.1        bit_4.5.0           
 [76] tidyselect_1.2.1     knitr_1.48           arrayhelpers_1.1-0   gridExtra_2.3        stats4_4.4.1        
 [81] xfun_0.47            bridgesampling_1.1-2 matrixStats_1.4.1    stringi_1.8.4        ggfun_0.1.6         
 [86] yaml_2.3.10          boot_1.3-30          evaluate_1.0.0       codetools_0.2-20     ggplotify_0.1.2     
 [91] cli_3.6.3            RcppParallel_5.1.9   Rdpack_2.6.1         xtable_1.8-4         munsell_0.5.1       
 [96] processx_3.8.4       coda_0.19-4.1        svUnit_1.0.6         parallel_4.4.1       rstantools_2.4.0    
[101] bayestestR_0.14.0    gap.datasets_0.0.6   bayesplot_1.11.1     Brobdingnag_1.2-9    lme4_1.1-35.5       
[106] mvtnorm_1.3-1        scales_1.3.0         insight_0.20.4       crayon_1.5.3         rlang_1.1.4         

Cite Packages

cite_packages()
  - Bürkner P (2017). "brms: An R Package for Bayesian Multilevel Models Using Stan." _Journal of Statistical Software_, *80*(1), 1-28. doi:10.18637/jss.v080.i01 <https://doi.org/10.18637/jss.v080.i01>. Bürkner P (2018). "Advanced Bayesian Multilevel Modeling with the R Package brms." _The R Journal_, *10*(1), 395-411. doi:10.32614/RJ-2018-017 <https://doi.org/10.32614/RJ-2018-017>. Bürkner P (2021). "Bayesian Item Response Modeling in R with brms and Stan." _Journal of Statistical Software_, *100*(5), 1-54. doi:10.18637/jss.v100.i05 <https://doi.org/10.18637/jss.v100.i05>.
  - Constantin A, Patil I (2021). "ggsignif: R Package for Displaying Significance Brackets for 'ggplot2'." _PsyArxiv_. doi:10.31234/osf.io/7awm6 <https://doi.org/10.31234/osf.io/7awm6>, <https://psyarxiv.com/7awm6>.
  - Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2024). _Rcpp: Seamless R and C++ Integration_. R package version 1.0.13, <https://CRAN.R-project.org/package=Rcpp>. Eddelbuettel D, François R (2011). "Rcpp: Seamless R and C++ Integration." _Journal of Statistical Software_, *40*(8), 1-18. doi:10.18637/jss.v040.i08 <https://doi.org/10.18637/jss.v040.i08>. Eddelbuettel D (2013). _Seamless R and C++ Integration with Rcpp_. Springer, New York. doi:10.1007/978-1-4614-6868-4 <https://doi.org/10.1007/978-1-4614-6868-4>, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). "Extending R with C++: A Brief Introduction to Rcpp." _The American Statistician_, *72*(1), 28-36. doi:10.1080/00031305.2017.1375990 <https://doi.org/10.1080/00031305.2017.1375990>.
  - Gabry J, Češnovar R, Johnson A, Bronder S (2024). _cmdstanr: R Interface to 'CmdStan'_. R package version 0.8.0.9000, commit a79cc5ef9355e9a9ec6a1539cab7fe69ec9b9917, <https://github.com/stan-dev/cmdstanr>.
  - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate." _Journal of Statistical Software_, *40*(3), 1-25. <https://www.jstatsoft.org/v40/i03/>.
  - Hartig F (2022). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. R package version 0.4.6, <https://CRAN.R-project.org/package=DHARMa>.
  - Hvitfeldt E (2021). _paletteer: Comprehensive Collection of Color Palettes_. R package version 1.3.0, <https://github.com/EmilHvitfeldt/paletteer>.
  - Kay M (2024). _tidybayes: Tidy Data and Geoms for Bayesian Models_. doi:10.5281/zenodo.1308151 <https://doi.org/10.5281/zenodo.1308151>, R package version 3.0.7, <http://mjskay.github.io/tidybayes/>.
  - Lenth R (2024). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. R package version 1.10.4, <https://CRAN.R-project.org/package=emmeans>.
  - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from Regression Models." _Journal of Open Source Software_, *3*(26), 772. doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
  - Makowski D, Lüdecke D, Patil I, Thériault R, Ben-Shachar M, Wiernik B (2023). "Automated Results Reporting as a Practical Tool to Improve Reproducibility and Methodological Best Practices Adoption." _CRAN_. <https://easystats.github.io/report/>.
  - Meredith M, Kruschke J (2022). _HDInterval: Highest (Posterior) Density Intervals_. R package version 0.2.4, <https://CRAN.R-project.org/package=HDInterval>.
  - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. R package version 3.2.1, <https://CRAN.R-project.org/package=tibble>.
  - Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, O'Hara R, Solymos P, Stevens M, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, De Caceres M, Durand S, Evangelista H, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill M, Lahti L, McGlinn D, Ouellette M, Ribeiro Cunha E, Smith T, Stier A, Ter Braak C, Weedon J (2024). _vegan: Community Ecology Package_. R package version 2.6-8, <https://CRAN.R-project.org/package=vegan>.
  - Pedersen T (2024). _patchwork: The Composer of Plots_. R package version 1.3.0, <https://CRAN.R-project.org/package=patchwork>.
  - R Core Team (2024). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.
  - Salazar G (2024). _EcolUtils: Utilities for community ecology analysis_. R package version 0.1, commit 2b57a8a8b76eba0852c56ce30f93d92c30320161, <https://github.com/GuillemSalazar/EcolUtils>.
  - Sarkar D (2008). _Lattice: Multivariate Data Visualization with R_. Springer, New York. ISBN 978-0-387-75968-5, <http://lmdvr.r-forge.r-project.org>.
  - Simpson G (2022). _permute: Functions for Generating Restricted Permutations of Data_. R package version 0.9-7, <https://CRAN.R-project.org/package=permute>.
  - Stan Development Team (2020). "StanHeaders: Headers for the R interface to Stan." R package version 2.21.0-6, <https://mc-stan.org/>.
  - Stan Development Team (2024). "RStan: the R interface to Stan." R package version 2.32.6, <https://mc-stan.org/>.
  - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_. Springer-Verlag New York. ISBN 978-3-319-24277-4, <https://ggplot2.tidyverse.org>.
  - Wickham H (2023). _forcats: Tools for Working with Categorical Variables (Factors)_. R package version 1.0.0, <https://CRAN.R-project.org/package=forcats>.
  - Wickham H (2023). _stringr: Simple, Consistent Wrappers for Common String Operations_. R package version 1.5.1, <https://CRAN.R-project.org/package=stringr>.
  - Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
  - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar of Data Manipulation_. R package version 1.1.4, <https://CRAN.R-project.org/package=dplyr>.
  - Wickham H, Henry L (2023). _purrr: Functional Programming Tools_. R package version 1.0.2, <https://CRAN.R-project.org/package=purrr>.
  - Wickham H, Hester J, Bryan J (2024). _readr: Read Rectangular Text Data_. R package version 2.1.5, <https://CRAN.R-project.org/package=readr>.
  - Wickham H, Vaughan D, Girlich M (2024). _tidyr: Tidy Messy Data_. R package version 1.3.1, <https://CRAN.R-project.org/package=tidyr>.
  - Wilke C (2024). _cowplot: Streamlined Plot Theme and Plot Annotations for 'ggplot2'_. R package version 1.1.3, <https://CRAN.R-project.org/package=cowplot>.
  - Yu G (2023). _ggimage: Use Image in 'ggplot2'_. R package version 0.3.3, <https://CRAN.R-project.org/package=ggimage>.
  - Zeileis A, Fisher JC, Hornik K, Ihaka R, McWhite CD, Murrell P, Stauffer R, Wilke CO (2020). "colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes." _Journal of Statistical Software_, *96*(1), 1-49. doi:10.18637/jss.v096.i01 <https://doi.org/10.18637/jss.v096.i01>. Zeileis A, Hornik K, Murrell P (2009). "Escaping RGBland: Selecting Colors for Statistical Graphics." _Computational Statistics & Data Analysis_, *53*(9), 3259-3270. doi:10.1016/j.csda.2008.11.033 <https://doi.org/10.1016/j.csda.2008.11.033>. Stauffer R, Mayr GJ, Dabernig M, Zeileis A (2009). "Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations." _Bulletin of the American Meteorological Society_, *96*(2), 203-216. doi:10.1175/BAMS-D-13-00155.1 <https://doi.org/10.1175/BAMS-D-13-00155.1>.